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We evaluate the free energy of the fluid and crystal phases for the ST2 potential [F.H. Stillinger 
and A. Rahman, J. Chem. Phys. 60, 1545 (1974)] with reaction field corrections for the long-range 
interactions. We estimate the phase coexistence boundaries in the temperature-pressure plane, 
as well as the gas-liquid critical point and gas-liquid coexistence conditions. Our study frames 
the location of the previously identified liquid-liquid critical point relative to the crystalline phase 
boundaries, and opens the way for exploring crystal nucleation in a model where the metastable 
liquid-liquid critical point is computationally accessible. 


I. INTRODUCTION 

The thermodynamic behavior of water at low temper¬ 
atures is unconventional. Several quantities, e.g. the 
isobaric density p, the isothermal compressibility Kt, 
and the constant-pressure specific heat Cp, are charac¬ 
terized by non-monotonic temperature or pressure de¬ 
pendence [l|. Over the past decades, the anomalous be¬ 
havior of these quantities has attracted the attention of 
numerous researchers. In 1992, a numerical investigation 
of the equation of state (EOS) s ugg ested the presence of 
a liquid-liquid (LL) critical point jj in the ST2 model Q, 
an interaction potential that describes water as a classi¬ 
cal, rigid, non-polarizable molecule. The presence of a LL 
critical point, located in the supercooled region, provides 
an elegant explanation of the thermodynamic anomalies 
that characterize liquid water and which become more 
pronounced close to such a critical point 0. 

The conceptual novelty of a one-component system 
with more than one liquid phase has stimulated the scien¬ 
tific communiW to deeply probe the physical origin of this 
phenomenon |5l-ll3l|. It is now clear that a LLCP, while 
common in tetrahedral network-forming liquids [l4l - ll9l| . 
can also be observed in complex one-component fluids 
when the (spherically symmetric) interaction potential 
generates two competing length scales (2nl - l^ . In the 
last few years the interest has shifted towards the inter¬ 
play between the liquid-liquid critical point and crystal 
nucleation (l^. [23 - It^ . Indeed, in experiments, crystal¬ 
lization has so far prevented direct observation of this 
phenomenon in a one-component bulk system. Only re¬ 
cently have computer simulations demonstrated the pos¬ 
sibility of generating a thermodynamically stable liquid- 
liquid critical point (as opposed to a metastable one) in 
models of network-forming liquids Bill. 

Accurate information on the phase coexistence bound¬ 
aries between disordered and ordered phases is relevant 
not only to establish the thermodynamic fields of sta¬ 
bility of the different phases, but also as a reference for 
estimating when the liquid becomes metastable. In turn, 
this has relevance for estimating when the barrier to crys¬ 


tallization becomes finite and how rapidly the barrier de¬ 
creases on supercooling [2^ . Except for one early report 
focussing on the liquid-ice Ih boundary [2^, none of the 
coexistence lines between the gas, liquid, and the many 
phases of crystalline ice have been accurately determined 
for the ST2 model. In this article we fill this gap and 
evaluate these coexistence boundaries by calculating the 
fluid chemical potential (via thermodynamic integration) 
and the crystal chemical potential (via the Frenkel-Ladd 
method , extended to molecules [HI). We test several 
crystals (ice lu, Ic, VI, VII, and VIII) and find that in 
the region of pressure where thermodynamic anomalies 
appear (e.g. near the lines of maxima of Cp and Kp) ice 
I;^ and Ic have the same free energy within our numeri¬ 
cal precision. Unexpectedly, we discover that for the ST2 
model, on increasing pressure, the stable phase is a dense 
tetragonal crystal with partial proton order. This struc¬ 
ture has a free energy about 0.4 kpT lower than ice VII, 
the structure obtained by interspersing two U lattices. 
(Here T is the temperature and kp is the Boltzmann 
constant.) We also evaluate the (metastable) line of co¬ 
existence for the recently reported ice 0 lattice [H, , a 

structure which could act (according to the Ostwald rule) 
as the intermediate phase in the process of nucleating 
the stable ice l^/c crystal from the fluid. For complete¬ 
ness, we determine the location of the gas-liquid critical 
point, which is found to be at Tf. = 558.0 ± 0.3K and 
Pc = 0.265 ± 0.005 g/cm3. 


II. MODEL AND SIMULATION METHODS 

We study, via Monte Carlo (MC) simulations, the 
original ST2 potential as defined by Rahman and Still¬ 
inger Q, with reaction field corrections to approximate 
the long-range contributions to the electrostatic interac¬ 
tions. ST2 models water as a rigid body with an oxygen 
atom at the center and four charges q = ±0.4e (where 
e is the electron charge), two positive and two nega¬ 
tive, in a tetrahedral geometry. The distances from the 
oxygen to the positive and negative charges are 0.1 and 


2 


0.08 nm respectively. The oxygen-oxygen interaction is 
modeled via a standard Lennard-Jones potential trun¬ 
cated at 2.5crLj, with a^j = 0.31 nm and e^j = 0.31694 
kJ/mol. The Lennard-Jones residual interactions are 
handled through standard long-range corrections, i.e. by 
assuming that the radial distribution function is unity 
beyond the cutoff. The charge-charge interactions are 
smoothly switched off both at small and large distances 
via a tapering function, as in the original model Q . Com¬ 
plete details of the simulation procedure are as described 
in Ref. Q. In the following, we use cr = 1 nm as unit of 
length. 



A. Thermodynamic integration: Fluid free energy 

To evaluate the fluid free energy we perform thermo¬ 
dynamic integration along a path of constant reference 
density p^ef for a modified pair potential, 

R = min(VsT 2 ,200 kJ/mol). (1) 

This potential coincides with the ST2 potential for all 
intermolecular distances and orientations where Vst 2 < 
200 kJ/mol, and is constant and equal to 200 kJ/mol 
otherwise. Note that in the temperature range where 
we investigate the phase behavior, molecules never ap¬ 
proach close enough to reach this limit. In this way, 
the divergence of the potential energy for configurations 
in which some intermolecular separations vanish (which 
would otherwise be probed at very high temperatures) is 
eliminated and the infinite temperature limit is properly 
approximated by an ideal gas of molecules at the same 
density. 

The fluid free energy (per particle) is calculated as 

rP 

Pref) = ^/ig(/3, Pref) + / p,ef)) d[3, (2) 

where /3 = l/fegT and /3/ig(;9,p) = log(p„cr^) - 1 is 
the ideal gas free energy and is the number density. 
Fig. [1] shows the average modified pair potential energy 
{V (/3, p)) and the interpolating (spline) continuous curve 
used to numerically evaluate the integral. The free energy 
at different densities along a constant-T path is evaluated 
via thermodynamic integration of the equation of state 

Pn) = p„,ref)+ r ^^^d\n{p'J, 

" Pn,Tei Pn 

(3) 

where P{pn) is the equation of state for the pressure P 
at fixed T. 


B. Crystal free energy 

To evaluate the free energy of a selected crys¬ 
talline structure we follow the methodology reviewed in 
Ref. [3l|. We define an Einstein crystal in which each 


FIG. 1: Average pair potential energy iV) vs. {RT)~^ at 
p = 1.0 g/cm® (with R the ideal gas constant). Symbols 
are MC data, and the line is the spline function used in the 
numerical integration. The inset shows the same data as a 
function of T and com par es them to previously published data 
for the ST2 potential . 


molecule interacts, in addition to the ST2 potential, with 
a Hamiltonian, composed of a translational (iJtrans) and 
a rotational (iJrot) part, that attaches each molecule to 
a reference position and orientation. For each particle 
we define two unit vectors: the (normalized) HH vector 
and dipole vector, named respectively a and b. The ref¬ 
erence configuration is defined by the reference position 
of the oxygen atom Tq and the reference position of a and 
b [3l|, [13 . In the following we indicate with r — rg the 
displacement of a particle located at r from its reference 
position, and with (j)a and (pb the angles between a and b 
and their reference values. More precisely. 


J^Einstein — Utrans T P-vot 


(4) 


with 


iJtrans = At(r-ro)V^" (5) 

and 


Prot 


— A 




( 6 ) 


Here At and A^ indicate the strength of the coupling to 
the reference configuration. Again following Ref. [ll| , the 
free energy (per particle) of a crystal structure /’'*, in the 
limit of large A^ and At is calculated as. 


Pr^ = Pfi+Ph + Ph+Ph + Ph + Ph (7) 










3 


where, indicating with N the number of molecules in the 
system, 


/3/i 




3(iV-l) 

2 


7V2 


1 

PnO-^ 


( 8 ) 


/3/2 

/3/3 

Pf4 

Ph 

Pk 


— In —1-1.5 ln(/3Ar) 

I iP^trans) X ^ 


{PHrot)xd^^P 

— ^VsT2 ^ 


ln(e“^^s'r2) 


Ar, At 


N 


( ln[1.5] (full proton-disordered crystal) 
( 0 (proton-ordered crystal). 


The symbols (Hrot)^ (Htrans);^ indicate the average 
values of iJrot and i?trans calculated from a MC simula¬ 
tion of particles interacting via the ST2 potential com¬ 
plemented by i?Einstein. The Symbol in¬ 

dicates the average value of e“^TsT 2 (where Vst 2 is the 
system ST2 potential energy) in a simulation in which the 
particles interact with each other via the ST2 potential 
and with the Einstein Hamiltonian with values A^. and A*. 
In all simulations carried out to perform the i nteg ration, 
the center of mass of the system is kept fixed 

Finally, Pk indicates the contribution of proton disor¬ 
der, evaluated according to Pauling’s estimate [33|. More 
recent calculations have essentially confirmed Pauling’s 
value 

Table U reports the values of Pfj for a few representa¬ 
tive cases. 


C. Grand canonical simulation: Gas-liquid phase 
coexistence 

To evaluate the gas-liquid coexistence and the loca¬ 
tion of the gas-liquid critical point, we perform grand- 
canonical MC simulations to evaluate at fixed T, volume 
V, and chemical potential p,, the probability p of observ¬ 
ing N particles in the simulated volume. To overcome the 
large free energy barriers separating the gas and liquid 
phases we implement the successive umbrella sampling 
(SUS) technique [1^. Since this method has been ap¬ 
plied previously to ST2 to estimate the liquid-liquid 
coexistence conditions, and has been documented in de¬ 
tail in these works, we refer the interested reader to the 
original literature. 


D. Proton position in the crystal structures 

To generate proton-disordered crystals, such as ice Ih/c 
and ice VII, one needs to assign protons to the oxygens. 



FIG. 2: Size effects associated with proton disorder in con¬ 
figurations with zero net dipole moment: Relation between 
the average pressure and the average energy of different 
proton-disordered configurations, for different number N of 
molecules. For ail cases, the structure is ice F at r=270 K 
and p = 0.8715 g/cnP. The inset shows the average pres¬ 
sure in different proton-disordered configurations for different 
number of molecules. 


located at the lattice positions, so as to satisfy the ice 
rules. To this end, we first calculate a list of all bonded 
oxygen neighbours (where four bonds connect to each 
oxygen atom) and then decorate the oxygen lattice by as¬ 
signing the proton for each bond to one of the two bonded 
atoms, iterating the following procedure: (i) Randomly 
select one oxygen with less than two hydrogens and one 
of the remaining undecorated bonds emanating from the 
selected oxygen, (ii) Randomly follow the path of un¬ 
decorated bonds until the path loops back to the original 
oxygen, (hi) Decorate all bonds of the selected path with 
one proton each, associating the protons to the oxygens 
encountered in the path. The procedure is iterated un¬ 
til all oxygens have two protons associated with them. 
Paths in which the initial and final oxygen atoms coin¬ 
cide only via periodic images produce a non-zero dipole 
moment and should be rejected if the net dipole moment 
of the cell is to vanish. 

To account for all possible proton realizations one 
needs to investigate large systems or average over several 
configurations. Indeed, we find that there is a signifi¬ 
cant correlation between the proton realization and the 
average potential energy E and average pressure P (at 
constant volume). Fig. [2] correlates P and E for each re¬ 
alization, while the inset shows P in different realizations 
for system sizes from N = 512 to iV = 21952 molecules. 
Only for 8000 or more particles is the variance between 
different realizations within a few MPa and a tenth of a 
kJ/mol, the tolerance required to allow for a precise de¬ 
termination of the thermodynamic variables entering into 
the free-energy calculation. Unless otherwise stated, we 
have analyzed configurations with 8000 or more particles 
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T(K) 

P (g/cm®) 

N 

Ar (kj/mol) 

At (kJ/mol) 

Pf 

Ph 

Ph 

Ph + PU 

Ph 

Ph 

Ice Ih 

270 

0.8715 

21952 

3- 

lO'^ 

3- 

10® 

-8.829 

19.440 

15.064 

-19.640 

-23.282 

-0.410 

Ice Ic 

270 

0.8715 

21952 

3- 

10"* 

3- 

10® 

-8.829 

19.440 

15.064 

-19.658 

-23.264 

-0.410 

Ice VI 

250 

1.27356 

8100 

3- 

10^ 

3- 

10® 

-10.772 

19.678 

15.305 

-20.349 

-25.996 

-0.410 

Ice VII 

270 

1.5804 

21296 

3- 

lO"* 

3- 

10® 

-8.230 

19.440 

15.064 

-19.318 

-23.006 

-0.410 

Ice vir 

270 

1.6250 

6912 

3- 

10^ 

3- 

10® 

-8.59 

19.440 

15.064 

-19.112 

-23.567 

-0.410 

Ice VIII 

270 

1.55645 

1152 

3- 

10^ 

3- 

10® 

-6.852 

19.4185 

15.064 

-19.061 

-22.274 

0 

Ice 0 

250 

0.8494 

29160 

3- 

10"* 

3- 

10® 

-10.399 

19.555 

15.180 

-19.983 

-24.741 

-0.410 

Fluid 

270 

1.002 

268 





-8.4411 

— 

— 

— 

— 

— 


TABLE I: Free energy /?/ of the fluid and crystal phases at selected points. The value for the residual entropy in Pf% was taken 
from Ref. [S^. The columns marked jSfi indicate the various contributions from the crystal free energy calculation. For the 
fluid, we used thermodynamic integration from an ideal gas at constant density, as explained in the text. 


for all proton-disordered crystals. 


III. RESULTS 
A. Gas-liquid coexistence 

Fig.Oshows the results of the SUS calculations. Panel 
(a) shows the probability p of finding N particles at fixed 
T and v at the coexistence chemical potential pc for dif¬ 
ferent T. Pc is evaluated by reweighting the histogram 
p{N) with respect to N, such that the area below the gas 
and the liquid peak is identical (0.5). At low T, the prob¬ 
ability minimum separating the two phases is more than 
50 orders of magnitude lower than the peak heights, high¬ 
lighting the need for a numerical technique (like SUS) 
that allows the observation of rare states. Close to the 
critical point [panel (b)], the probability of exploring in¬ 
termediate densities between the gas and the liquid be¬ 
comes significant and p{N) [or p{p)] assumes the char¬ 
acteristic shape typical of all systems belonging to the 
same universality class. Panel (c) compares p{N + sE), 
where E is the potential energy of the configuration and 
s is the so-called mixing field parameter |41| . with the 
theoretical expression for the magnetization in the Ising 
model. To reinforce the identification of the critical point 
with the Ising universality class, the inset shows the finite 
size scaling of the critical T (defined as the T, for each 
size, at which the fluctuations in N + sE are best fitted 
with the Ising form) as a function of 
with 9 = 0.54 and v = 0.630 [1^, H^. The extrapola¬ 
tion to L —>■ oo suggests that the gas-liquid critical point 
for the reaction field ST2 model is Tc = 558.0 ± 0.3K 
and Pc = 0.265 ± 0.005 g/cm®. Finally, panel (d) shows 
the gas-liquid coexistence in the p — T plane. A clear 
nose appears around T = 300 K, signaling the onset of 
the network of hydrogen bonds (HB). Indeed, strong di¬ 
rectional interactions (such as the HB), impose a strong 
coupling between density and energy. The formation of a 
fully bonded tetrahedral network (the expected thermo¬ 
dynamically stable state at low T) requires a well-defined 
minimum local density, which for the present model is ap¬ 


proximately p = 0.8 g/cm®. Hence, at low T, the density 
of the network coexisting with the gas must approach this 
value. For completeness, the inset in panel (d) reports 
the value of Ppc along the coexistence line. 


B. Fluid-crystal coexistence 

We have investigated the stability of crystal phases 
that may coexist with the fluid at low T. In particular, 
we have determined the free energies of ices U, U, VI, 
VB, and VBI, as well as the recently proposed metastable 
ice 0 structure [s^- Note that with the exception of ice 
VBI, all these phases have disordered hydrogen bonding. 
Examples of our thermodynamic integration results are 
reported in Fig. |4j where we plot the reduced chemical 
potential /3p = (Sf + pP/pn of different phases at two se¬ 
lected T. For each pressure interval, the lowest chemical 
potential phase is the thermodynamically stable one. In¬ 
tersections of different curves locate coexistence points, 
either stable or metastable. We then interpolate the fluid 
and crystal free energies based on the equation of state 
to draw the coexistence lines in the phase diagram. 

The complete phase diagram is reported in Fig. (5] At 
low T and low P, the most stable crystal structure is the 
ice I lattice. From our simulations, the cubic (Ic) and 
hexagonal (U) ice structures have the same free energy 
within our numerical accuracy. At positive pressures, the 
liquid phase coexisting with ice I is always denser than 
ice, and as a result, the melting temperature of ice I de¬ 
creases with increasing P. At negative P (near P = —80 
MPa), the ice I and liquid phases coexist at the same den¬ 
sity, and the melting temperature reaches a maximum. 
We note that we have confirmed the ice I/j/c melting 
temperature calculated via thermodynamic integration 
at two separate pressures using direct coexistence simu¬ 
lations, and find good agreement. We note that for the 
ST2-Ewald model, the melting temperature of U at the 
single pressure of 260 MPa was estimated to be around 
274 K, consistent with the present estimate [^ . 

At high pressure, the main candidate structures are 
the proton-ordered ice VBI structure, and the proton- 
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FIG. 3: Gas-liquid coexistence for the ST2 model, (a) Distribution of the density fluctuations in gas-liquid coexisting states for 
different T for system size L = 2 nm. (b) Same data as in (a) but for T close to the critical temperature Tc and for L = 6 nm. 
(c) Gomparison of the fluctuation in the order parameter x (a linear combination of N and E) with the theoretical expression 
for the three-dimension Ising model. The inset shows the finite size scaling of the critical temperature, (d) Resulting gas-liquid 
phase diagram in the T — p plane. The inset shows the values of the chemical potential along the coexistence. 


disordered ice VII structure. Both structures consist of 
two interpenetrating Ic lattices (somewhat distorted in 
the case of proton-ordering), where the oxygen positions 
form a BCC lattice. According to our free-energy cal¬ 
culations, the disordered ice VII is the more stable one 
in the region where coexistence with the fluid might oc¬ 
cur. However, when trying to confirm the accuracy of our 
predicted liquid-ice VII coexistences using direct coexis¬ 
tence simulations, we observed crystal growth at temper¬ 
atures significantly above the melting temperature pre¬ 
dicted from free energy calculations. The newly grown 
parts of the crystal still display the BCC topology of the 
oxygen atoms, but the crystal shrinks by a few percent 
in the direction perpendicular to the growth direction, 
leading to a slight distortion of the lattice, that we refer 
to in the following as ice VII*. As this distortion does 
not occur in fully disordered ice VII, we attribute the 
unexpectedly high stability of the ice VII* lattice to the 
emergence of partial proton ordering, which decreases the 
crystal free energy. To confirm this, we created a fully re¬ 
grown ice VII* configuration by alternately melting and 
regrowing the two halves of an ice VII configuration in an 


elongated simulation box. When measuring the proton- 
proton and dipole-dipole correlation functions for both 
the original ice VII structure and the regrown ice VII*, 
we see only minor changes in the proton-proton correla¬ 
tion function in the region 3 A< r < 4 A [see Fig. Ha)]. 
In contrast, the dipole-dipole correlation function [see 
Fig. Hb)] shows significant additional signal which al¬ 
though weak, extends up to long spatial scales. Using 
the Frenkel-Ladd method, we calculate the free energy 
of this configuration (assuming full proton disorder), and 
find that it is indeed lower than that of the original crys¬ 
tal by Ri 0.4 ksT per particle, confirming that the lower 
melting temperature observed in our direct coexistence 
simulations can be attributed to the (slight) change in 
crystal structure. The difference in free energy mainly 
results from the lower potential energy of the regrown 
crystal. We note here that partial proton ordering would 
reduce the contribution of the residual entropy to the 
free energy of the crystal, causing us to underestimate 
the ice VII* free energy. On the other hand, the pres¬ 
ence of defects in the system is expected to cause an 
overestimate in the crystal free energy. It is thus not a 
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FIG. 4: Reduced chemical potential /3/i as a function of pres¬ 
sure P for competing phases at r=270 K and T=300 K. At 
each pressure, the phase with the lowest chemical potential 
is the stable one. Crossings indicate (metastable or stable) 
phase transitions. 




priori obvious that this free energy can be used to predict 
coexistences. Nonetheless, comparing the melting tem¬ 
perature predicted from the free energy and equations 
of state of the regrown crystal with the melting temper¬ 
ature taken from the direct coexistence simulations, we 
find good agreement (T « 320K ± 5K at P = 250 MPa). 
Calculating the rest of the coexistence lines for this crys¬ 
tal using thermodynamic integration, we observe that ice 
VII* has a significantly larger stability region than the 
original ice VII (see Fig. [5]). 

We note that neither ice VI nor ice 0 are ever the most 
thermodynamically stable phase in the investigated re¬ 
gion. As it may be relevant in future nucleation studies, 
we include the metastable coexistence line of the liquid 
with ice 0 in the phase diagram (Fig. [S]). 


FIG. 5: (a) Pressure-temperature phase diagram for the ST2 
model with reaction field. In this P-T region, the stable 
phases are the ice VII*, cubic or hexagonal ice (ice R/c), 
gas, and liquid. The thick lines indicate phase boundaries be¬ 
tween stable phases, while the thinner lines denote metastable 
phase transitions. The full circle indicates the location of the 
liquid-liquid critical point and its error bars, as estimated 
in Ref. (b) Temperature-density representation of the 

same phase diagram. White regions indicate two-phase coex¬ 
istence. Filled regions indicate areas of one-phase stability for 
the different phases. Solid lines indicate coexistence densities. 
Note that the one-phase stability field of ice l^/c is centred 
on the optimal crystal density and is very narrow, compara¬ 
ble in width to the coexistence lines. Horizontal dotted lines 
correspond to triple points. We use the same colour-coding as 
in (a). The blue dots indicate the coexisting densities of the 
low and high density liquids, from Ref. [^. The red square 
is the estimated location of the LL critical point. 
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r (nm) 


FIG. 6: (a) Radial distribution function and (b) dipole dipole 
correlation function for the fully proton-disordered ice VII 
structure and the regrown ice VII* crystal, at p = 1.673 
g/cm®, in the inherent structure. Note the slight extra cor¬ 
relation in gunir) and gdip-dip(r) in the region 2 A< r < 
4 A. 


IV. CONCLUSIONS 

Recently, the ST2 potential has been at the centre of 
renewed interest in connection to the debate on the origin 
of the liquid-liquid critical point P, . This model 

exhibits known deficiencies in accurately modelling water 
properties, e.g. it overemphasizes the tetrahedrality of 
the liquid structure, thus shifting all water anomalies to 


higher temperatures. Despite these deficiencies, the ST2 
model plays a key role as a prototype system in many 
studies related to the presence of a liquid-liquid critical 
point. We report here fundamental properties of the ST2 
model, by evaluating the location of the gas-liquid critical 
point and the gas-liquid coexistence curve, as well as the 
coexistence lines between the liquid and several crystal 
structures, allowing us to map out the phase diagram of 
the ST2 model in the low-temperature regime. We find a 
stable ice I phase at low pressure and temperature, with 
both the hexagonal and cubic stackings approximately 
equal in free energy. Differently from real water, the 
high-pressure phase behavior of the model is dominated 
by a new crystal whose growth is templated by the ice VII 
interface. This ice VII* tetragonal crystal is composed of 
a lattice in which the oxygens have the same topology as 
ice VII but in which the protons are not completely ran¬ 
domly distributed. We have not been able to identify a 
small unit cell for this new crystal, but inspection of the 
HH radial distribution function indicates minute but ob¬ 
servable differences in the region around 3.3 A, accompa¬ 
nied by weak but long ranged correlations in the dipole- 
dipole correlation function. This structure, despite the 
small partial proton order, has a significant lower poten¬ 
tial energy than VII (approximately 1.2 kJ/mol). As a 
result, ice VII* is significantly more stable than the fully 
proton-disordered ice VII phase at all pressures and it 
dominates the high-pressure phase behavior of the model. 
The liquid-liquid critical point for this model lies, ac¬ 
cording to the most recent estimates, inside the region of 
stability of the ice VII* crystal phase and is metastable 
with respect to ice R or R as well as to ice VII. Our 
results provide a starting point for the study of nucle- 
ation in the ST2 model, as well as for the exploration 
of modifications to the model that could make the 
liquid-liquid critical point more accessible [l8l| . 


V. ACKNOWLEDGMENTS 

We thank L. Filion, A. Geiger, A. Rehtanz, and I. 
Saika-Voivod for useful discussions. We are honored to 
dedicate this study to Prof. Jean-Pierre Hansen, from 
whom we have learned liquid state theory. 


[1] P. G. Debenedetti and H. E. Stanley, Phys. Today 56, 
40 (2003). 

[2] P. H. Poole, F. Sciortino, U. Essmann, and H. E. Stanley, 
Nature 360, 324 (1992). 

[3] F. H. Stillinger and A. Rahman, J. Chem. Phys. 60, 1545 
(1974). 

[4] L. Xu, P. Kumar, S. V. Buldyrev, S.-H. Chen, P. H. 
Poole, F. Sciortino, and H. E. Stanley, Proc. Natl. Acad. 
Sci. U.S.A. 102, 16558 (2005). 


[5] O. Mishima and H. E. Stanley, Nature 396, 329 (1998). 

[6] A. K. Soper and M. A. Ricci, Phys. Rev. Lett. 84, 2881 

( 2000 ). 

[7] Y. Katayama, T. Mizutani, W. Utsumi, O. Shimomura, 
M. Yamakata, and K.-i. Funakoshi, Nature 403, 170 
( 2000 ). 

[8] R. Kurita and H. Tanaka, Science 306, 845 (2004). 

[9] A. Taschin, P. Bartolini, R. Eramo, R. Righini, and 
R. Torre, Nat. Commun. 4 (2013). 





































[10] G. Pallares, M. E. M. Azouzi, M. A. Gonzalez, J. L. 
Aragones, J. L. Abascal, G. Valeriani, and F. Gaupin, 
Proc. Natl. Acad. Sci. U.S.A. Ill, 7936 (2014). 

[11] K. Amann-Winkel, C. Gainaru, P. H. Handle, M. Seidl, 
H. Nelson, R. Bohmer, and T. Loerting, Proc. Natl. 
Acad. Sci. U.S.A. 110, 17720 (2013). 

[12] M. E. M. Azouzi, C. Ramboz, J.-F. Lenain, and 
F. Gaupin, Nature Phys. 9, 38 (2013). 

[13] J. A. Sellberg, C. Huang, T. McQueen, N. Loh, 
H. Laksmono, D. Schlesinger, R. Sierra, D. Nordlund, 
C. Hampton, D. Starodub, et al.. Nature 510 , 381 (2014). 

[14] L Saika-Voivod, F. Sciortino, and P. H. Poole, Phys. Rev. 
E 63, 011202 (2001). 

[15] V. V. Vasisht, S. Saw, and S. Sastry, Nat. Phys. 7, 549 

( 2011 ). 

[16] C. W. Hsu, J. Largo, F. Sciortino, and F. W. Starr, Proc. 
Natl. Acad. Sci. U.S.A. 105 , 13711 (2008). 

[17] J. L. Abascal and C. Vega, J. Ghem. Phys. 133, 234502 

( 2010 ). 

[18] F. Smallenburg, L. Filion, and F. Sciortino, Nature Phys. 
10, 653 (2014). 

[19] F. W. Starr and F. Sciortino, Soft Matter 10, 9413 
(2014). 

[20] E. Jagla, J. Ghem. Phys. Ill, 8980 (1999). 

[21] G. Franzese, G. Malescio, A. Skibinsky, S. V. Buldyrev, 
and H. E. Stanley, Nature 409, 692 (2001). 

[22] L. Xu, S. V. Buldyrev, G. A. Angell, and H. E. Stanley, 
Phys. Rev. E 74, 031108 (2006). 

[23] P. Gallo and F. Sciortino, Phys. Rev. Lett. 109, 177801 

( 2012 ). 

[24] D. T. Limmer and D. Ghandler, J. Ghem. Phys. 135 , 
134503 (2011). 

[25] J. G. Palmer, F. Martelli, Y. Liu, R. Car, A. Z. Pana- 
giotopoulos, and P. G. Debenedetti, Nature 510 , 385 
(2014). 

[26] R. S. Singh and B. Bagchi, J. Ghem. Phys. 140 , 164503 
(2014). 

[27] C. R. C. Buhariwalla, R. K. Bowles, 1. Saika-Voivod, 
F. Sciortino, and P. H. Poole, preprint, arxiv:1501.03115 
(2015). 

[28] F. Romano, E. Sanz, and F. Sciortino, J. Ghem. Phys. 
134 , 174502 (pages 8) (2011). 

[29] A preliminary study of the liquid-ice R coexistence line 


for the ST2 model was carried out by A. Rehtanz, A. 
Geiger, and P.H. Poole; see A. Rehtanz, Ph.D. thesis, 
Dortmund University (Logos Verlag, Berlin, 2000). 

[30] D. Frenkel and A. J. C. Ladd, J. Ghem. Phys 81, 3188 
(1984). 

[31] C. Vega, E. Sanz, J. Abascal, and E. Noya, J. Phys.: 
Condens. Matter 20 , 153101 (2008). 

[32] J. Russo, F. Romano, and H. Tanaka, Nature materials 
13 , 733 (2014). 

[33] D. Quigley, D. Alf, and B. Slater, J. Ghem. Phys. 141, 
161102 (2014). 

[34] P. H. Poole, S. R. Becker, F. Sciortino, and F. W. Starr, 
The Journal of Physical Chemistry B 115 , 14176 (2011). 

[35] E. G. Noya, M. M. Conde, and C. Vega, J. Ghem. Phys. 
129 , 104704 (2008). 

[36] B. Smith and D. Frenkel, Understanding molecular sim¬ 
ulations (Academic Press, New York, 1996). 

[37] L. Pauling, The Nature of the Chemical Bond (Cornell 
University Press, Ithaca, New York, 1945). 

[38] B. A. Berg, C. Muguruma, and Y. Okamoto, Phys. Rev. 
B 75 , 092202 (2007). 

[39] P. Virnau and M. Muller, J. Ghem. Phys. 120, 10925 
(2004). 

[40] F. Sciortino, 1. Saika-Voivod, and P. H. Poole, Phys. 
Ghem. Ghem. Phys. 13 , 19759 (2011). 

[41] N. B. Wilding, Phys. Rev. E 52 , 602 (1995). 

[42] A. M. Ferrenberg and D. Landau, Phys. Rev. B 44, 5081 
(1991). 

[43] A. Pelissetto and E. Vicari, Phys. Rep. 368 , 549 (2002). 

[44] M. J. Cuthbertson and P. H. Poole, Phys. Rev. Lett. 106, 
115706 (2011). 

[45] G. Franzese, G. Malescio, A. Skibinsky, S. V. Buldyrev, 
and H. Stanley, Nature 409, 692 (2001). 

[46] F. Sciortino, E. La Nave, and P. Tartaglia, Phys. Rev. 
Lett. 91 , 155701 (2003). 

[47] D. Fuentevilla and M. Anisimov, Phys. Rev. Lett. 97 , 
195702 (2006). 

[48] V. Holten, C. E. Bertrand, M. A. Anisimov, and J. V. 
Sengers, J. Ghem. Phys. 136, 094507 (2012). 

[49] V. Holten, J. C. Palmer, P. H. Poole, P. G. Debenedetti, 
and M. A. Anisimov, J. Ghem. Phys. 140, 104502 (2014). 

[50] F. Smallenburg and F. Sciortino, to be published (2015). 


